Chapter 6 - Compaction and its inherent length scale¶
The compaction-press problem¶
%matplotlib inline
import numpy as np
import scipy.sparse as sps
import scipy.sparse.linalg as spla
from scipy.optimize import brentq
import matplotlib.pyplot as plt
from matplotlib import animation
import warnings
warnings.filterwarnings('ignore')
plt.rcParams["animation.html"] = "jshtml"
plt.ioff();
The solution to the Filter Press problem was given as
\[
\begin{equation}
\label{eq:filterpress-cmprate}
\cmp(z) = -\frac{W_0}{\cmplength_0}\,\e^{z/\cmplength_0} \quad \text{ for } \quad z\le0.
\end{equation}
\]
def compaction_rate(z):
"""
Solution to the Filter Press problem
:param z:
:return: the normalized solution of the Filter Press problem
"""
return -np.exp(z)
The compact rate function, equation \(\eqref{eq:filterpress-cmprate}\), is plotted in Figure below.
z = np.linspace(-8, 0.0, 1000)
C = compaction_rate(z)
fig, ax = plt.subplots()
fig.set_size_inches(18.0, 4.5)
ax.plot([0., 0], [-2., 2], '--k', linewidth=1)
ax.plot([-8.0, 1.0], [0., 0], ':k', linewidth=1)
ax.plot(z, C, '-k', linewidth=3)
ax.set_xlabel('$z/\delta_0$', fontsize=20)
ax.set_ylabel('$\mathcal{C}\,\delta_0/W_0$', fontsize=20)
ax.text(-7.0, -0.2, '$W_0 \Rightarrow$', fontsize=20)
ax.set_xlim(-8.0, 0.3)
ax.set_ylim(-1.1, 0.1)
ax.tick_params(axis='both', which='major', labelsize=13)
plt.show()
The permeability-step problem¶
The porosity is given by the piece-wise constant function,
\[\begin{split}
\begin{equation}
\label{eq:permstep-porosity}
\por(z) =
\begin{cases}
f_p\por_0 & \text{for $z > 0$},\\
f_m\por_0 & \text{for $z \le 0$},
\end{cases}
\end{equation}
\end{split}\]
where \(f_i\) (\(i=p,m\)) are constants that multiply the reference porosity, chosen such that \(f_i-\vert 1\vert\ll 1\).
The solution of the permeability-step problem was given as
\[\begin{split}
\begin{equation}
\label{eq:permstep-solution-full}
\cmp(z) = \frac{\por_0w_0}{\cmplength_0}\tilde{f}
\begin{cases}
\exp\left(-z/\cmplength_p\right) & \text{for $z>0$},\\
\exp\left(z/\cmplength_m\right) & \text{for $z\le 0$}.
\end{cases}
\end{equation}
\end{split}\]
def cmprate(fp, fm, f0, n, z):
"""
solution of the permeability-step problem
:param fp:
:param fm:
:param f0:
:param n:
:param z:
:return:
"""
cmp = (np.power(fm, n) * (1.0 - fm*f0)
- np.power(fp, n) * (1.0 - fp*f0))/(np.power(fm, 0.5*n) + np.power(fp, 0.5*n))
return np.asarray([cmp*np.exp(-z_* np.power(fp, -0.5*n)) if z_ > 0.0
else cmp*np.exp(z_ * np.power(fm, 0.5*n)) for z_ in z])
The one-dimensional segregation flux \(q \equiv \phi(w-W)\) is given by
\[\begin{split}
\begin{equation}
\label{eq:permstep-darcy-flux-solution}
q = \por_0w_0
\begin{cases}
\left[f_p^\permexp(1-f_p\por_0) + \tilde{f}f_p^{\permexp/2}\exp(-z/\cmplength_p)\right] & z>0,\\[2mm]
\left[f_m^\permexp(1-f_m\por_0) - \tilde{f}f_m^{\permexp/2}\exp(z/\cmplength_m)\right] & z\le0.
\end{cases}
\end{equation}
\end{split}\]
def segflux(fp, fm, f0, n, z):
"""
one-dimensional segmentation flux
:param fp:
:param fm:
:param f0:
:param n:
:param z:
:return:
"""
cmp = (np.power(fm, n) * (1.0 - fm*f0)
- np.power(fp, n) * (1.0 - fp*f0))/(np.power(fm, 0.5*n) + np.power(fp, 0.5*n))
return np.asarray([np.power(fp, n)*(1.0-fp*f0) + cmp*np.power(fp, 0.5*n)*np.exp(-z_* np.power(fp, -0.5*n))
if z_ > 0.0
else np.power(fm, n)*(1.0 - fm*f0) - cmp* np.power(fm, 0.5*n)*np.exp(z_* np.power(fm, -0.5*n))
for z_ in z])
Case 1¶
The porosity increases with z: \(f_m=0.85\), \(f_p=1.15\).
fm, fp = 0.85, 1.15
z = np.linspace(-4.0, 4.0, 1000)
f = np.asarray([fm if z_<0.0 else fp for z_ in z])
fig, ax = plt.subplots()
fig.set_size_inches(4.5, 9.0)
ax.plot(f, z, '-k', linewidth=2)
ax.plot(np.ones(1000), z, ':k', linewidth=1)
ax.set_ylabel(r'$z/\delta_0$', fontsize=20)
ax.set_xlabel(r'$\phi/\phi_0$', fontsize=20)
ax.tick_params(axis='both', which='major', labelsize=13)
plt.show()
The figure below shows that the scaled compaction rate and segregation flux for the porosity increase:
zmax = 4
n = 3
f0 = 0.01
z = np.linspace(-4.0, 4.0, 1000)
C = cmprate(fp, fm, f0, n, z)
fwmW = segflux(fp, fm, f0, n, z)
fig, ax = plt.subplots()
fig.set_size_inches(4.5, 9.0)
ax.plot(C, z, '-k', linewidth=2, label='$\mathcal{C}(z)\delta_0/(\phi_0w_0)$')
ax.plot(fwmW, z,'--k', linewidth=2, label='$q(z)/(\phi_0w_0)$')
ax.set_ylabel(r'$z/\delta_0$', fontsize=20)
ax.legend(fontsize=15, loc='upper center')
ax.tick_params(axis='both', which='major', labelsize=13)
plt.show()
Case 2¶
The porosity decreases with z: \(f_m=1.15\), \(f_p=0.85\)
fm, fp = 1.15, 0.85
f = np.asarray([fm if z_<0.0 else fp for z_ in z])
fig, ax = plt.subplots()
fig.set_size_inches(4.5, 9.0)
ax.plot(f, z, '-k', linewidth=2)
ax.plot(np.ones(1000), z, ':k', linewidth=1)
ax.set_xticks((0.85, 1.0, 1.15))
ax.set_ylabel(r'$z/\delta_0$', fontsize=20)
ax.set_xlabel(r'$\phi/\phi_0$', fontsize=20)
ax.tick_params(axis='both', which='major', labelsize=13)
plt.show()
The figure below shows that the scaled compaction rate and segregation flux for the porosity decrease:
zmax = 4
n = 3
f0 = 0.01
z = np.linspace(-4.0, 4.0, 1000)
C = cmprate(fp, fm, f0, n, z)
fwmW = segflux(fp, fm, f0, n, z)
fig, ax = plt.subplots()
fig.set_size_inches(4.5, 9.0)
ax.plot(C, z, '-k', linewidth=2, label='$\mathcal{C}(z)\delta_0/(\phi_0w_0)$')
ax.plot(fwmW, z,'--k', linewidth=2, label='$q(z)/(\phi_0w_0)$')
ax.tick_params(axis='both', which='major', labelsize=13)
ax.legend(fontsize=20, loc='lower center')
plt.show()
Propagation of small porosity disturbances¶
The phase (\(c_p\)) and group (\(c_g\)) velocities are given as
\[\begin{split}
\begin{align}
\label{eq:linwave-phase-vel}
c_p &= \frac{\permexp}{k^2 + 1},\\
c_g &= c_p - \frac{2\permexp k^2}{(k^2 + 1)^2}.
\end{align}
\end{split}\]
The phase (\(c_p\)) and group (\(c_g\)) velocities are plotted below as a function of the wavelength \(\lambda=2\pi/k\).
lambdas = np.logspace(-2.0, 3.0, 1000)
k = 2*np.pi/lambdas
n = 2.0
cp_2 = n/(k ** 2 + 1)
cg_2 = cp_2 - 2.0 * n * (k**2) / ((k**2 + 1.0) ** 2)
n = 3.0
cp_3 = n/(k**2 + 1)
cg_3 = cp_3 - 2.0 * n * (k**2) / ((k**2 + 1.0) ** 2)
fig, ax = plt.subplots()
fig.set_size_inches(9.0, 9.0)
ax.semilogx(lambdas, cp_2, '--k', linewidth=3, label='$c_p, n=2$')
ax.semilogx(lambdas, cg_2, '--k', linewidth=1, label='$c_g, n=2$')
ax.semilogx(lambdas, cp_3, '-k', linewidth=3, label='$c_p, n=3$')
ax.semilogx(lambdas, cg_3, '-k', linewidth=1, label='$c_g, n=3$')
ax.plot([10**(-5), 10**5], [0.0, 0.0], '-', color=[0.75, 0.75, 0.75])
ax.plot([2.0*np.pi, 2.0*np.pi],[-10.0, 10.0],'-', color=[0.75, 0.75, 0.75])
ax.set_xlabel(r'$\lambda$', fontsize=20)
ax.set_ylabel(r'$c$', fontsize=20)
ax.set_ylim(-0.5, 3.1)
ax.set_xlim(10**(-2), 10**3)
ax.tick_params(axis='both', which='major', labelsize=13)
plt.legend(fontsize=20)
plt.show()
Magmatic solitary waves¶
The solitary wave speed \(v\) is computed as
\[\begin{split}
\begin{equation*}
\solwavespeed(\solwaveamp) =
\begin{cases}
(\solwaveamp-1)^2/\left(\solwaveamp\ln\solwaveamp - \solwaveamp+1\right)& \text{for }\permexp=2,\\
2\solwaveamp + 1 & \text{for }\permexp=3.
\end{cases}
\end{equation*}
\end{split}\]
The dimensionless solitary wave speed \(\upsilon\) computed with the equation above as a function of wave amplitude \(\Lambda^*\) relative to the background porosity is plotted below.
lambdas = np.linspace(1.001, 4.0, 1000)
c2 = (lambdas - 1.0)**2 / (lambdas * np.log(lambdas) - lambdas + 1.0)
c3 = 2.0 * lambdas + 1.0
fig, ax = plt.subplots()
fig.set_size_inches(9.0, 9.0)
ax.plot(lambdas, c2, '--k', linewidth=2, label='$n=2$')
ax.plot(lambdas, c3, '-k', linewidth=2, label='$n=3$')
ax.set_xlim(1.0, 4.0)
ax.set_ylim(0.0, 7.0)
ax.set_xlabel(r'$\Lambda^*$', fontsize=20)
ax.set_ylabel(r'$\upsilon$', fontsize=20)
ax.tick_params(axis='both', which='major', labelsize=13)
ax.legend(fontsize=20)
plt.show()
The normalized porosity profile is obtained through numerical inversion of
\[
\begin{equation}
\label{eq:magmons-implict-shape}
Z(\solwave) = \left(\solwaveamp+\tfrac{1}{2}\right)^{\tfrac{1}{2}}\left[
\frac{1}{(\solwaveamp-1)^{\tfrac{1}{2}}}\ln\left(
\frac{(\solwaveamp-1)^{\tfrac{1}{2}}-(\solwaveamp-\solwave)^{\tfrac{1}{2}}}
{(\solwaveamp-1)^{\tfrac{1}{2}}+(\solwaveamp-\solwave)^{\tfrac{1}{2}}}
\right)-2\left(\solwaveamp-\solwave\right)^{\tfrac{1}{2}}\right]
\end{equation}
\]
which is implemented in Python as
def porosity(f, z, A):
"""
normalized porosity profile
:param f:
:param z:
:param A:
:return:
"""
sqrtAf = np.sqrt(A-f)
sqrtAm1 = np.sqrt(A - 1.)
return z + np.sqrt(A+0.5)*(-2.*sqrtAf + np.log((sqrtAm1-sqrtAf)/(sqrtAm1+sqrtAf))/sqrtAm1)
Profiles of normalised porosity perturbation for solitary waves of various amplitude are plotted below. We set \(n=3\) in all cases.
zmax = 25.0
zs = np.linspace(0.0, zmax, 100)
zm = 0.5*(zs[1:]+zs[0:-1])
AsLS = {2.: '-.', 4.: '-', 6.:'--', 8.: ':'}
phi = {}
for A, ls in AsLS.items():
phi[A] = np.asarray([brentq(lambda f: porosity(f, z, A), 1.000000001, A) for z in zs])
fig, ax = plt.subplots()
fig.set_size_inches(18.0, 9.0)
fig.set_facecolor('w')
for A, ls in AsLS.items():
ax.plot(np.concatenate((-zs[::-1], zs), axis=0), np.concatenate((phi[A][::-1], phi[A]), axis=0),
'k', label='$\Lambda^* = '+str(int(A))+'$', linestyle=ls, linewidth=2)
ax.set_xlim(-zmax, zmax)
ax.set_ylim(0.0, 9.0)
ax.grid()
ax.set_xlabel(r'$Z$, compaction lengths', fontsize=20)
ax.set_ylabel(r'$\Lambda$', fontsize=20)
ax.tick_params(axis='both', which='major', labelsize=13)
ax.legend(fontsize=20)
plt.show()
Profiles of compaction rate for solitary waves of various amplitude. \(n=3\) in all cases. The gravity vector \(\gravity\) points to the left, as indicated in the figure below.
zmax = 25.0
zs = np.linspace(0.0, zmax, 100)
zm = 0.5*(zs[1:]+zs[0:-1])
AsLS = {2: '-.', 4: '-', 6:'--', 8: ':'}
C = {}
for A, ls in AsLS.items():
C[A] = -(2.0*A + 1.0)*(phi[8][1:]-phi[8][0:-1])/(zs[1]-zs[0]) # phi is defined in the previous cell
fig, ax = plt.subplots()
fig.set_size_inches(18.0, 9.0)
for A, ls in AsLS.items():
plt.plot(np.concatenate((-zm[::-1], zm), axis=0), np.concatenate((-C[A][::-1], C[A]), axis=0),
'k', linestyle=ls, linewidth=2)
ax.set_xlim(-zmax, zmax)
ax.set_ylim(-10.0, 10.0)
ax.grid()
ax.set_xlabel(r'Z, compaction lengths', fontsize=20)
ax.set_ylabel(r'$\lambda$', fontSize=20)
ax.tick_params(axis='both', which='major', labelsize=13)
plt.show()
Solitary-wave trains¶
The equations below admit a nonlinear solitary wave solution:
\[\begin{split}
\begin{align}
-\left(\npor^\permexp\cmp_z\right)_z + \cmp &= -\left(\npor^\permexp\right)_z, \\
\npor_t &= \cmp.
\end{align}
\end{split}\]
which is plotted below:
def get_compaction_rate_dirichlet(phi, n, dz, phi0):
n_ = len(phi)
perm = np.sqrt((phi[0:-1] ** n) * (phi[1:] ** n))
# rhs
b = np.zeros(n_, dtype=float)
b[1:-1] = phi0 * dz * (perm[1:] - perm[0:-1])
# matrix
offsets = np.array([0, -1, 1])
data = np.zeros(3*n_).reshape(3, n_)
data[0, 0] = data[0, -1] = 1
data[0, 1:-1] = -(perm[0:-1] + perm[1:] + dz * dz) # diagonal
data[1, 0:-2] = perm[0:-1] # sub-diagonal
data[2, 2:] = perm[1:] # sup-diagonal
mtx = sps.dia_matrix((data, offsets), shape=(n_, n_))
mtx = mtx.tocsr()
# solution of linear system
Cmp = spla.dsolve.spsolve(mtx, b)
return Cmp
def solitary_wave_update_porosity(PhiOld, n, phi0, dz, dt):
Cmp = get_compaction_rate_dirichlet(PhiOld, n, dz, phi0)
PhiNew = PhiOld + dt * Cmp / phi0
Cmp = 0.5 * (Cmp + get_compaction_rate_dirichlet(PhiNew, n, dz, phi0))
PhiNew = PhiOld + dt * Cmp / phi0
return PhiNew
fig, ax = plt.subplots(figsize=(4.5, 9.0))
ln, = plt.plot([], [], 'k')
phi0 = 0.05 # background porosity
A = 1.5 # amplitude of step
zmax = 150. # total size of domain
z0 = zmax / 5. # location of step
zw = 10. # width of step
n = 3. # permeability exponent
Nz = 1000 # number of grid points
cfl = 1. # courant limit on time-step
tmax = 70. # maximum time
# initial condition
z = np.linspace(0.0, zmax, Nz)
f = 1. - (A - 1.) * (1 + np.tanh((z - z0) / zw)) / 2.
# derived parameters
dz = z[1] - z[0]
V = (1. - phi0 ** n) / (1. - phi0)
dt = cfl * dz / V
Nt = int(np.ceil(tmax / dt))
def init():
ax.set_xlim(0.4, 1.5)
ax.set_ylim(0., zmax)
ax.tick_params(axis='both', which='major', labelsize=13)
return ln,
def update(frame):
global f
f = solitary_wave_update_porosity(f, n, phi0, dz, dt)
ln.set_data(f, z)
return ln,
animation.FuncAnimation(fig, update, frames=np.linspace(0, tmax, Nt), init_func=init, blit=True)